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A computationally simple way to accommodate 'basins' of trapping states in standard kinetic 
Monte Carlo simulations is presented. By assuming the system is effectively equilibrated in the 
basin, the residence time (time spent in the basin before escape) and the probabilities for transition 
to states outside the basin may be calculated. This is demonstrated for point defect diffusion over 
a periodic grid of sites containing a complex basin. 



The kinetic Monte Carlo (kMC) method is used to 
evolve atomistic systems dynamically from state to state 
over timescales much longer than can be achieved in 
molecular dynamics simulations [l], 0] ■ The method uti- 
lizes a catalog of state-to-state transition rates obtained 
from atomistic (dynamic or static) calculations, to de- 
termine probabilistically a sequence of states (and their 
residence times) that closely resembles the actual sys- 
tem dynamics. The computational efficiency of the kMC 
method is due to the neglect of details: the system is 
simply moved from one distinct state to another, and 
the time clock is advanced accordingly. However, it may 
be that the set of transition rates is such as to equili- 
brate the system in a subset of mutually accessible states, 
from which escape is a very rare event. This situation of 
course reduces the efficiency of the method greatly. Here 
we present a simple means to accommodate such equi- 
librating basins in the standard kMC approach for the 
case of defect diffusion in solids or on surfaces. In fact 
the basin is regarded as just another accessible defect site 
with a characteristic residence time. This is only possi- 
ble when the defect is considered to have equilibrated in 
the basin (that is, all sites in the basin have been visited 
many times), so its entry and exit points are uncorre- 
lated. This treatment of equilibrating basins will be par- 
ticularly useful for kMC simulations of defect diffusion in 
nanocrystalline materials, where the diffusion coefficients 
for the defect in the grain boundaries and the crystalline 
grains may differ by many orders of magnitude and 
of radiation damage in solids, where microstructure evo- 
lution (driven by defect diffusion) over very long time 
scales is of interest. 

In a kMC simulation of defect diffusion, the defect 
moves over a regular or irregular grid of sites (represent- 
ing the potential wells that can accommodate the defect) 
according to probabilistic rules. The diffusion coefficient 
D is then obtained in the usual way: D = {x'^)/{2dt), 
where x is the defect displacement over the time t, and 
d is the dimension of the space. Typically, the resi- 
dence times associated with moves from visited sites are 
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summed until the required time interval t is completed. 
But for purposes of this derivation, it is necessary to also 
regard t as the sum of the accrued residence times at the 
visited sites. This is because those accrued times, for sites 
in the basin, are proportional to the equilibrium defect 
concentrations there. Also for purposes of the derivation, 
it is convenient to use the term 'periphery site' for those 
sites in the basin from which the defect can move out of 
the basin. With this terminology set, we obtain expres- 
sions for the probability pi of escape from the basin via 
the particular periphery site i, and for the residence time 
^basin associated with a visit by the defect to the basin. 

Consider a defect at periphery site i. With each 
move from site i, the accrued residence time for 
that site increases by an (average) amount xt = 

(Efc(i) h^b{i) + J2q{i) h^q(t) ) > where the two sums are 
over all transition rates from site i to accessible sites b{i) 
within the basin and to accessible sites q{i) outside the 
basin, respectively. The probability that it will escape 

the basin on that move is = thus 

on average one of every visits by the defect to site i 
will result in an escape. In that event, the residence time 

tf = r^e^^ = {j2q{i)ki->q{i)^ , on average, has accrued 
to site i. It is noteworthy that tf is a function only of the 
rates fcj^q(j) out of the basin. 

Of course, the basin may contain many periphery and 
interior sites. As the defect moves within the basin, it 
produces an increasingly accurate set {tk/{tk)} of rela- 
tive residence times, where the sites k are in the basin 
(periphery and interior) and the average (indicated by 
the angle brackets) is taken over all sites in the basin. 
In fact the elements tk/{tk) approach the values Ck/{ck), 
where Cfc is the equilibrium defect concentration at site 
k that is routinely obtained in molecular dynamics and 
statics calculations (cfc is an exponential function of the 
defect formation energy at site k). During the time 
T = X^fe^fei the number of visits by the defect to site i 
is U/Ti. Since the probability that a particular visit will 
not lead to an escape from the basin is (1 — e^), the a pri- 
ori probability that the defect does not escape from the 
basin via site i during the time interval T is (1 — e^)*'/'^'. 
Then the a priori probability that the defect does es- 
cape the basin via site i during the time interval T is 
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1 - (1 - ej)*'/""' « {U/Ti)ei for e, < 1. This equals U/t<r 
when the definition = Tie~^ is used. Thus the proba- 
bility Pi that the defect escapes the basin from periphery 
site i rather than from another periphery site is given by 



(1) 



where the sum is over all periphery sites j . Substituting 
into Eq. ([T|) the expression for tf gives 
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As evident from this last equation, the escape from pe- 
riphery site i out of the basin would be to site q' with 

probability pi^g' = h^q' (Eg(i)fcj^g(o) ' where site q' 
is one of the set {<?(«)}. That is, a defect trapped in the 
basin will escape via periphery site i to site q' (outside 
the basin) with probability Pi^q' — PiPi^q' . 

The long-term, average behavior of the defect is thus 
reproduced by the standard kMC method, with the addi- 
tion that if the defect enters the basin, on its next move 
it escapes the basin from a periphery site chosen in ac- 
cordance with the probability distribution implied by Eq. 

The residence time tbasin associated with this move 
is given by the relation 
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where the first sum is over all periphery sites j, and the 
second sum is over all basin interior sites k. The sec- 
ond sum accounts for the time the defect spends at in- 
terior sites, which in the case of a particular interior site 
k equals the ratio tk/tj of time spent by the defect at 
site k to time spent at an arbitrarily chosen periphery 
site J, multiplied by the average time f^pj spent at site 
J during a visit by the defect to the basin (note that 
the ratio tjPj/tj is identical for all periphery sites j). 
The simplest example demonstrating Eq. ([3]) is that of 
a basin comprised of n identical periphery sites (that is, 
the probability pj of escaping the basin via a particular 
periphery site is 1/n, and all equal the 'lifetime' t^) 
and no interior sites. Clearly the average residence time 
in the basin per visit by the defect is t^, so the average 
residence time in each periphery site per visit to the basin 
must be i'^/n, which equals fjpj as expected. 

By use of Eq. ^ for pj , Eq. ^ may be rewritten as 
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where now the sum is over all (periphery and interior) 
basin sites k. As discussed above, the ratio tk/tj may be 
replaced by Ck /cj . Thus the equilibrating basin is accom- 
modated by addition of the set {Pi~tq(i) \ of probabilities 



for moves out of the basin, and the residence time ^basin, 
to the kMC catalog of transition rates. 

It may be noted that the derivation of Eq. ^ relies 
on the use of the average value (called r above) for the 
time that accrues to a basin site with each visit by the 
defect prior to escape. In conventional kMC simulations, 
the time may alternatively be advanced by an amount 
At taken randomly from the exponential distribution 
T^"'" exp(— At/r); that is, by the amount At = r[— Inz] 
where z is chosen randomly from the interval (0, 1]. Thus 
it is possible in the latter case to calculate the higher mo- 
ments of the escape time from the basin as well as the 
average time ibasin- Of course, the method developed 
here for handling deep basins in kMC simulations pre- 
supposes that calculation of an accurate distribution of 
basin escape times (whether desired or not) is not com- 
putationally feasible. In this event, it is recommended 
(for consistency) that average values r, rather than vari- 
able values At, be used to accrue time to sites outside 
the basin. This should not affect the average value {x^) 
obtained for a specified diffusion time t, that is needed 
to calculate the defect diffusion coefficient D. 

This method of handling a set of connected states may 
be contrasted with that of Novotny Q , who applies the 
finite Markov chain formalism The basin sites are 
therefore transient states, and the sites to which the de- 
fect moves out of the basin are absorbing states. All 
transition probabilities connecting transient states, and 
connecting transient states with absorbing states, are el- 
ements in the Markov transition matrix M. Then the 
formalism gives, for the defect in a specified initial tran- 
sient state, (1) the mean number of times in each of the 
transient states before absorption, and (2) the probabil- 
ities for absorption in each of the absorbing states. (See 
Ref. ^ for a detailed example of how to use finite Markov 
chain theory to model stochastic physical systems.) The 
correlation between the entrance and exit points at the 
basin periphery is thus preserved at the expense of con- 
siderable mathematical and computational complication 
(e.g., a different matrix M is needed for each of the pos- 
sible initial states). That virtue is minor when the defect 
is essentially equilibrated in the basin before its escape, 
and in any event may be negated by the various sources of 
error (e.g., inaccurate transition rates) and the stochastic 
nature of the simulation. It should be emphasized that 
the Markov approach requires that all transition rates be- 
tween basin sites be available, while the present approach 
can alternatively use equilibrium defect concentrations. 

Before applying the method to sample systems with 
complex basins, it is interesting to consider a very simple, 
one-dimensional system that can be solved analytically. 
This is a linear arrangement of four sites, labeled (in 
order) 1 through 4, where the transition rates fe^a and 
k^^2 are much faster than the rates k2^i and k^^n. Thus 
a defect will 'flicker' between sites 2 and 3 many times 
before escaping to site 1 or 4 Q. The average behavior 
of the defect in this system is easily calculated by use of 
the Markov formalism when sites 1 and 4 are regarded as 
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absorbing states. In the event that the defect is initially 
at site 2, the analytic calculation produces the row vector 

P2^1 + P2^3P3->4: 

where Pi^j is the probability for the defect at site i to 
move to site j [so, for example, P2^i — fc2-»i 7(^2^1 + 
^2^3)]; the elements f3i and /?2 are the probabilities for 
absorption at site 1 and site 4, respectively; and the ele- 
ments (3^ and P4 are the mean number of times at sites 
2 and 3, respectively, before absorption. The expressions 
for Pi and P2 have been obtained previously by Mason 
et al. Q , by accounting for all possible numbers of flick- 
ers prior to escape from sites 2 and 3; for example, the 
probability that a defect initially at site 2 will escape to 
site 1 is Y,'^=o{P2->3P3^2Tp2^i ==^2^1/(1-^2^3^3^2), 
which equals (3i. 

In the event that the defect is initially at site 3, the 
corresponding calculation produces the row vector 

P' = — (P3^2P2^1 P3^4 P3^2 l) ■ 

P2~>1 +P2->3P3^i 

Then the 'averaged' results are given by the row vector 
/? = X2/3 + X3/5', where X2 and X3 are relative concen- 
trations at sites 2 and 3 that satisfy X2 + X3 = 1 and 
detailed balance, X2^2^3 = X3^3^2- Note that this av- 
eraging removes any memory of the 'initial' defect site 
(that is, whether the defect entered from site 1 or from 
site 4) . The averaged vector is 



P2^1 +P2^3P3^4 
{llP3^2P2~,l 72P2^3P3^4 71^3^2 72^2^3) 

where 71 = 1 -I- /c3__4(fc2^3 + ^3^2)""^ and 72 = 1 + 
fc2_>i(fc2— ►s -l- fc3_>2)~^. This may be compared with the 
equivalent row vector B constructed from the stochastic 
quantities derived above for an equilibrated basin: 

B = (P2 P3 

1 

P3^2P2^1 + P2^3P3^4 
{P3^2P2^1 P2^3P3^A P3->2 P2— s) 

which very closely resembles (3 when P2-+3 S> P2^i and 

A more complex basin is represented in Fig. 1. This 
system is a periodically repeated (in both dimensions) 
10 X 10 regular network of nodes (defect-accessible sites) 
connected by bonds (diffusion paths), where the 'equi- 
librating basin' is the subset of 34 nodes connected by 
the 40 thick bonds. Given the transition rates associated 
with each bond, it is a straightforward matter to obtain 
the defect diffusion coefficient by a kMC simulation. 

Table I presents the diffusion coefficients calculated 
by the standard method ('Exact') and by the 'basin' 
method ('Approx.'), and an estimate of the relative 



FIG. 1: Representation of a system of trapping and non- 
trapping sites. Those sites (nodes) connected by the thick 
bonds comprise the equilibrating basin in which the defect 
may be trapped for very long periods of time. 



computation time needed in each case, for three dif- 
ferent sets of transition rates. The first set (row 1) 
has ki^j — 10exp[— (/ii — ^j)] for the thick bonds and 
ki^j = exp[— (/Zi— /J,j)] for the thin bonds, where the {^i} 
are chemical potentials assigned to the nodes with values 
taken randomly from the interval [0, 1]. The second set 
(row 2) is similar to the first set, but with the difference 
that the for nodes belonging to the basin are taken 
from the interval [3,4], so that the defect will segregate 
to the basin. The third set (row 3) is similar to the first 
set, but with the rates ki^j for the thick bonds having 
the prefactor 1000 (instead of 10). With these transi- 
tion rates, detailed balance is satisfied: Ciki^j = Cjkj^i. 
The set {q} is needed to calculate the probabilities {pi} 
and the residence time tbasin, and furthermore provides a 
nice check on the calculations (namely, the accrued res- 
idence time ti at node i should be proportional to Ci). 
The values for the diffusion coefficient D are believed to 
be accurate to ±1 in the last digit. In the last column, 
the 'speed-up factor' (due to use of the basin method) 
refers to the computational time needed to accomplish a 
given defect diffusion time t, not to the computational 
time needed to achieve a particular accuracy. 

The large difference in D values in the first row of 
Table I shows that the basin method does a poor job 
when the defect cannot equilibrate before escaping; that 
is, when there is a significant spatial correlation between 
the entry and exit points (in this case due to the small 
diffusivity contrast between regions, which does not suf- 
ficiently confine the defect to the basin). Otherwise, the 
diffusion coefficients obtained by assuming the defect to 
equilibrate in the basin are seen to be very comparable to 
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Transition rates 


Exact D 


Approx. D 


Speed-up factor 


^(th^ck) ^ 10cxp[-(/ii - ^j)] 








k'^J^f^ = exp[-(//i - i^j)] 


1.373 


1.799 


1.2 


iJ.i e [0, 1] 








^(thick) ^ iOexp[-(/i, - //j)] 
A:'!!;^' =exp[-(/Xi -Mj)] 

^(basin) ^ p^^j 


0.0285 


0.0287 


3.1 


^(non-ba.in) ^ jq^ 








^(thick) ^ iooOexp[-(^, -^,)] 








kf^J^ =exp[-(/xi -/xj)] 


1.79 


1.799 


81.0 


Mi € [0, 1] 









TABLE I: Comparison of diffusion coefficients calculated by the kinetic Monte Carlo method. The transition rates k are for 
the paths represented by thick and thin bonds in Fig. 1; the /i are chemical potentials associated with the nodes that, for the 
purposes of this work, ensure that detailed balance is obeyed. The 'Exact' and 'Approx.' D are diffusion coefficients calculated 
in the standard manner, and with the set of trapping states treated as an equilibrating basin, respectively. The 'Speed-up 
factor' shows the computational advantage of the latter approach. 



the 'exact' values, while costing (potentially) orders-of- 
magnitude less computer time. Furthermore, the accrued 
residence times at the nodes (both inside and outside the 
basin) are in every case extremely close to their exact 
values (proportional to the {cj}). 

The results in Table I give a general indication of the 
utility of the basin method. In particular, the method 
is accurate when the defect is essentially equilibrated in 
the basin. The extent to which this is the case may be 
determined by a conventional kMC simulation (where the 
basin method is not used): the set {tk/{tk)} of relative 
residence times for sites in the basin, obtained for a single 
visit by the defect to the basin, is compared with the set 
{cfe/(cfe)}. The two sets are more or less identical for a 
defect that is more or less equilibrated in the basin. 

In general the basin method gives an upper bound for 
the actual diffusion coefficient. This is due to its ne- 
glect of any spatial correlation between the entry and 
exit points at the basin periphery: the distance be- 
tween these points is, on average, less when they are 
spatially correlated than when they are not. In either 
case the time spent in the basin per visit has average 
value <basin (calculated according to the analytic expres- 
sion above), so a higher value for the diffusion coefficient 
is obtained in the latter case. [That the average time 
spent in the basin per visit is tbasin m both cases is evi- 



dent from the fact that a kMC simulation will produce 
a set {tm/{tm)} ~ {cm/{cm)} (where now all sites m 
in the system — those outside the basin as well as those 
inside — are included), whether the basin method is incor- 
porated in the kMC code or not.] A comparison of rows 
1 and 3 in Table I illustrates this point. The two systems 
with different sets of transition rates nonetheless possess 
(by design) identical sets {c„i}, {pj}, and {kj^q(j)}, and 
identical basin residence time tbasin: this is the reason 
the two systems produce the same 'Approx.' value for 
the diffusion coefficient (1.799). But the defect in the 
first system (row 1) is not well equilibrated in the basin, 
causing an 'Approx.' value for D that is too high in that 
case. 

As a final comment, it should be emphasized that this 
approach to accommodating such trapping basins (cre- 
ated by, for example, segregation or orders-of-magnitude 
differences in transition rates as considered in Table I) 
in kMC simulations gives increasingly accurate results as 
the degree of confinement increases, which is precisely 
the situation where kMC simulations are, in the absence 
of this approach, increasingly inefficient and inaccurate. 
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